home *** CD-ROM | disk | FTP | other *** search
/ CD Ware Multimedia 1995 May / cd Ware (Juegos) Epimundo.iso / DOS / C / MATV.ZIP / MATV.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-11  |  19.0 KB  |  856 lines
  1.  
  2. // Matv - A Simple Matrix Class
  3.  
  4. // Version: 1.0
  5. // Author: Mark Von Tress, Ph.D.
  6. // Date: 10/11/92
  7.  
  8. // Copyright(c) Mark Von Tress 1992
  9.  
  10.  
  11. // DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  12. // WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  13. // TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  14. // ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  15. // FROM USE OF THIS PROGRAM.
  16.  
  17.  
  18. #include "matv.h"
  19.  
  20. //////////////////////////////// vector
  21.  
  22. ///// error handler
  23. FLOAT Vector::Nrerror(const char *errormsg)
  24. {
  25.    FLOAT x = 2;
  26.    cout << "\nMATRIX ERROR: " << errormsg << "\nexiting to system\n";
  27.    cerr << "\nMATRIX ERROR: " << errormsg << "\nexiting to system\n";
  28.    x = 2;
  29.    exit(1);
  30.    return x;
  31. }
  32.  
  33. ///// Matrix integrity
  34. void Vector::Garbage(const char *errormsg)
  35. {
  36.    int errorint = 0;
  37.    if (signature != SIGNATURE) errorint++;
  38.    if (n <= 0) errorint++;
  39.    if (c[-1] != signature) errorint++;
  40.    if (errorint) {
  41.       cout << "Operating on a Corrupted matrix: " << "\n";
  42.       Nrerror(errormsg);
  43.    }
  44. }
  45.  
  46. ///// index
  47. FLOAT &Vector::operator () (int i)
  48. {
  49.    return c[i];
  50. }
  51.  
  52. ///// output
  53. ostream &operator << (ostream &s, Vector &v)
  54. {
  55. #ifndef NO_CHECKING
  56.    v.Garbage("Vector <<: v is garbage");
  57. #endif
  58.    s << "CanDelete: " << (int) v.CanDelete << " CanAlias: "
  59.      << (int) v.CanAlias << " c:\t" << v.c << "\n";
  60.    for (int i = 0; i < v.n; i++) {
  61.       s << v(i) << " ";
  62.       if (!(i % 7) && i > 0) s << "\n";
  63.    }
  64.    s << "\n";
  65.    return s;
  66. }
  67.  
  68. ///// memory allocator
  69. FLOAT *Vector::SetupVectors(int nn)
  70. {
  71.    static unsigned long myconst =
  72.           (65536 / ((unsigned long) sizeof (FLOAT))) - 2;
  73.  
  74.    if (nn <= 0) Nrerror("SetupVectors: n < 0");
  75.    if (nn > myconst)
  76.       Nrerror("SetupVectors: more than 64K in vector(n>16382)");
  77.  
  78.    n = nn;
  79.    c = new FLOAT[n + 1];
  80.    if (!c) Nrerror("SetupVectors: allocation error");
  81.    c[0] = signature = SIGNATURE;
  82.    c++;
  83.    c[0] = 1;
  84.    return c;
  85. }
  86.  
  87. ///// memory deletion
  88. void Vector::PurgeVectors(void)
  89. {
  90.    signature++;
  91.    if (*-- c == SIGNATURE) {
  92.       c[0] = 0;
  93.       delete[] c;
  94.    }
  95.    else
  96.       Vector::Nrerror("PurgeVectors: attempting to delete a deleted vector");
  97. }
  98.  
  99. ////////// constructors
  100. ///// default constructor
  101. Vector::Vector(void) : n(1), CanDelete(TTRUE), CanAlias(TTRUE),
  102.         signature(SIGNATURE)
  103. {
  104.    c = SetupVectors(1);
  105. }
  106.  
  107. ///// allocate an n vector
  108. Vector::Vector(int nn) : n(nn), CanDelete(TTRUE), CanAlias(TTRUE),
  109.         signature(SIGNATURE)
  110. {
  111.    c = SetupVectors(nn);
  112. }
  113.  
  114. ///// copy an array into a vector
  115. Vector::Vector(int nn, FLOAT *x) : n(nn), CanDelete(TTRUE),
  116.         CanAlias(TTRUE), signature(SIGNATURE)
  117. {
  118.    c = SetupVectors(nn);
  119.    FLOAT *cc = c;
  120.    FLOAT *acc = x;
  121.    int i = n;
  122.    while (i--) *cc++ = *acc++;
  123. }
  124.  
  125. ///// copy constructor
  126. Vector::Vector(Vector &a) : n(a.n), CanDelete(TTRUE), CanAlias(TTRUE),
  127.         signature(SIGNATURE)
  128. {
  129.    if (a.CanDelete == FFALSE && a.CanAlias) {
  130.       // a is not responsible for deleting c and a can have an alias
  131.       c = a.c;
  132.    }
  133.    else {
  134.       c = SetupVectors(a.n);
  135.       FLOAT *cc = c;
  136.       FLOAT *acc = a.c;
  137.       int i = n;
  138.       while (i--) *cc++ = *acc++;
  139.    }
  140. }
  141.  
  142. ///// destructor
  143. Vector::~Vector()
  144. {
  145.    n = 0;
  146.    if (CanDelete == TTRUE)
  147.       PurgeVectors();
  148. }
  149.  
  150. ////// assignment
  151. Vector &Vector::operator = (Vector &a)
  152. {
  153. #ifndef NO_CHECKING
  154.    a.Garbage("Vector Assignment: a is garbage");
  155.    Garbage("Vector Assignment: *this is garbage");
  156. #endif
  157.    if (this == &a) return *this;
  158.    // deal with aliased storage
  159.    if (c == a.c && a.CanDelete == FFALSE) {
  160.       // *this is aliased with a, and a is not
  161.       // responsible for deleting c.
  162.       CanDelete = TTRUE;
  163.       n = a.n;
  164.       return *this;
  165.    }
  166.    if (c == a.c && a.CanDelete == TTRUE) {
  167.       // *this is aliased with a, and a is
  168.       // responsible for deleting c.
  169.       c = SetupVectors(a.n);
  170.       CanDelete = TTRUE;
  171.    }
  172.    // *this is no longer aliased with a
  173.    if (n != a.n) {
  174.       PurgeVectors();
  175.       c = SetupVectors(a.n);
  176.    }
  177.    FLOAT *cc = c;
  178.    FLOAT *acc = a.c;
  179.    int i = n;
  180.    while (i--) *cc++ = *acc++;
  181.    return *this;
  182. }
  183.  
  184. //////////////////////// matrix class
  185.  
  186. ///// indexing
  187. FLOAT &Matrix::operator () (int i, int j) {
  188. #ifndef NO_CHECKING
  189.    if (i < 1 || i > r || j < 1 || j > c)
  190.       Nrerror("operator(): index out of range");
  191. #endif
  192.    return Vector::operator () ((i - 1) *c + (j - 1));
  193. }
  194.  
  195. ///// assignment
  196. Matrix &Matrix::operator = (Matrix &a)
  197. {
  198. #ifndef NO_CHECKING
  199.    a.Garbage("Matrix Assignment: a is garbage");
  200.    Garbage("Matrix Assignment: *this is garbage");
  201. #endif
  202.    r = a.r;
  203.    c = a.c;
  204.    this->Vector::operator = (a);
  205.    return *this;
  206. }
  207.  
  208. ///// output
  209. ostream &operator << (ostream &s, Matrix &m)
  210. {
  211. #ifndef NO_CHECKING
  212.    m.Garbage("<<: m is garbage");
  213. #endif
  214.    s << "\n";
  215.    for (int i = 1; i <= m.r; i++) {
  216.       for (int j = 1; j <= m.c; j++)
  217.          s << m(i, j) << " ";
  218.       s << "\n";
  219.    }
  220.    s << "\n";
  221.    return s;
  222. }
  223.  
  224. void Matrix::show(char *str)
  225. {
  226. #ifndef NO_CHECKING
  227.    Garbage("show: *this is garbage");
  228. #endif
  229.    cout << str << "\n";
  230.    cout << *this;
  231. }
  232.  
  233. //////////// matrix operators
  234.  
  235. ///// addition
  236. Matrix operator + (Matrix &a, Matrix &b)
  237. {
  238. #ifndef NO_CHECKING
  239.    a.Garbage("Matrix +: a is garbage");
  240.    b.Garbage("Matrix +: b is garbage");
  241. #endif
  242.    if (a.r != b.r || a.c != b.c)
  243.       a.Nrerror("matrices do not conform in addition");
  244.    
  245.    Matrix c(a.r, a.c);
  246.    for (int i = 1; i <= c.r; i++)
  247.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) + b(i, j);
  248.    c.release();
  249.    return c;
  250. }
  251.  
  252. Matrix operator + (Matrix &a, FLOAT b)
  253. {
  254. #ifndef NO_CHECKING
  255.    a.Garbage("Matrix +: a is garbage");
  256. #endif
  257.    
  258.    Matrix c(a.r, a.c);
  259.    for (int i = 1; i <= c.r; i++)
  260.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) + b;
  261.    c.release();
  262.    return c;
  263. }
  264.  
  265. Matrix operator + (FLOAT b, Matrix &a)
  266. {
  267. #ifndef NO_CHECKING
  268.    a.Garbage("Matrix +: a is garbage");
  269. #endif
  270.  
  271.    Matrix c(a.r, a.c);
  272.    for (int i = 1; i <= c.r; i++)
  273.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) + b;
  274.    c.release();
  275.    return c;
  276. }
  277.  
  278. Matrix &Matrix::operator += (Matrix &a)
  279. {
  280. #ifndef NO_CHECKING
  281.    a.Garbage("Matrix +=: a is garbage");
  282.    Garbage("Matrix +=: *this is garbage");
  283. #endif
  284.    if (a.r != r || a.c != c)
  285.       a.Nrerror("matrices do not conform in addition");
  286.  
  287.    for (int i = 1; i <= r; i++)
  288.       for (int j = 1; j <= c; j++) (*this)(i, j) +=  a(i,j);
  289.  
  290.    return *this;
  291. }
  292.  
  293. Matrix &Matrix::operator += (FLOAT b)
  294. {
  295. #ifndef NO_CHECKING
  296.    Garbage("Matrix +=: *this is garbage");
  297. #endif
  298.  
  299.    for (int i = 1; i <= r; i++)
  300.       for (int j = 1; j <= c; j++) (*this)(i, j) +=  b;
  301.  
  302.    return *this;
  303. }
  304.  
  305. ///// subtraction
  306. Matrix operator - (Matrix &a, Matrix &b)
  307. {
  308. #ifndef NO_CHECKING
  309.    a.Garbage("Matrix -: a is garbage");
  310.    b.Garbage("Matrix -: b is garbage");
  311. #endif
  312.    if (a.r != b.r || a.c != b.c)
  313.       a.Nrerror("matrices do not conform in subtraction");
  314.  
  315.    Matrix c(a.r, a.c);
  316.    for (int i = 1; i <= c.r; i++)
  317.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) - b(i, j);
  318.    c.release();
  319.    return c;
  320. }
  321.  
  322. Matrix operator - (Matrix &a, FLOAT b)
  323. {
  324. #ifndef NO_CHECKING
  325.    a.Garbage("Matrix -: a is garbage");
  326. #endif
  327.  
  328.    Matrix c(a.r, a.c);
  329.    for (int i = 1; i <= c.r; i++)
  330.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) - b;
  331.    c.release();
  332.    return c;
  333. }
  334.  
  335. Matrix operator - (FLOAT b, Matrix &a)
  336. {
  337. #ifndef NO_CHECKING
  338.    a.Garbage("Matrix -: a is garbage");
  339. #endif
  340.  
  341.    Matrix c(a.r, a.c);
  342.    for (int i = 1; i <= c.r; i++)
  343.       for (int j = 1; j <= c.c; j++) c(i, j) = b - a(i, j);
  344.    c.release();
  345.    return c;
  346. }
  347.  
  348. Matrix operator - (Matrix &a)
  349. {
  350. #ifndef NO_CHECKING
  351.    a.Garbage("Matrix Unary -: a is garbage");
  352. #endif
  353.    Matrix c(a.r, a.c);
  354.    for (int i = 1; i <= c.r; i++)
  355.       for (int j = 1; j <= c.c; j++) c(i, j) = -a(i, j);
  356.    c.release();
  357.    return c;
  358. }
  359.  
  360. Matrix &Matrix::operator -= (Matrix &a)
  361. {
  362. #ifndef NO_CHECKING
  363.    a.Garbage("Matrix -=: a is garbage");
  364.    Garbage("Matrix -=: *this is garbage");
  365. #endif
  366.    if (a.r != r || a.c != c)
  367.       a.Nrerror("matrices do not conform in subtraction");
  368.  
  369.    for (int i = 1; i <= r; i++)
  370.       for (int j = 1; j <= c; j++) (*this)(i, j) -=  a(i,j);
  371.  
  372.    return *this;
  373. }
  374.  
  375. Matrix &Matrix::operator -= (FLOAT b)
  376. {
  377. #ifndef NO_CHECKING
  378.    Garbage("Matrix -=: *this is garbage");
  379. #endif
  380.  
  381.    for (int i = 1; i <= r; i++)
  382.       for (int j = 1; j <= c; j++) (*this)(i, j) -=  b;
  383.  
  384.    return *this;
  385. }
  386.  
  387. ///// multiplication
  388. Matrix operator *(Matrix &a, Matrix &b)
  389. {
  390. #ifndef NO_CHECKING
  391.    a.Garbage("Matrix *: a is garbage");
  392.    b.Garbage("Matrix *: b is garbage");
  393. #endif
  394.    if (a.c != b.r)
  395.       a.Nrerror("matrices do not conform in multiplication");
  396.  
  397.    Matrix c(a.r, b.c);
  398.    LFLOAT sum;
  399.    for (int i = 1; i <= a.r; i++)
  400.       for (int j = 1; j <= b.c; j++) {
  401.          sum = 0;
  402.          for (int k = 1; k <= a.c; k++)
  403.             sum += (LFLOAT) ((LFLOAT) a(i, k)) *((LFLOAT) b(k, j));
  404.          c(i, j) =(FLOAT) sum;
  405.    }
  406.    c.release();
  407.    return c;
  408. }
  409.  
  410. Matrix operator *(Matrix &a, FLOAT b)
  411. {
  412. #ifndef NO_CHECKING
  413.    a.Garbage("Matrix *: a is garbage");
  414. #endif
  415.  
  416.    Matrix c(a.r, a.c);
  417.    for (int i = 1; i <= c.r; i++)
  418.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) *b;
  419.    c.release();
  420.    return c;
  421. }
  422.  
  423. Matrix operator *(FLOAT b, Matrix &a)
  424. {
  425. #ifndef NO_CHECKING
  426.    a.Garbage("Matrix *: a is garbage");
  427. #endif
  428.  
  429.    Matrix c(a.r, a.c);
  430.    for (int i = 1; i <= c.r; i++)
  431.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) *b;
  432.    c.release();
  433.    return c;
  434. }
  435.  
  436. Matrix &Matrix::operator *= (Matrix &a)
  437. {
  438. #ifndef NO_CHECKING
  439.    a.Garbage("Matrix *=: a is garbage");
  440.    Garbage("Matrix *=: *this is garbage");
  441. #endif
  442.    if (a.r != c)
  443.       a.Nrerror("matrices do not conform in multiplication");
  444.  
  445.    LFLOAT sum = 0;
  446.    Matrix cc( r, a.c );
  447.    for (int i = 1; i <= cc.r; i++)
  448.       for (int j = 1; j <= cc.c; j++){
  449.           sum = 0;
  450.           for( int k=1; k<=c; k++)
  451.               sum += (LFLOAT) ((LFLOAT) (*this)(i, k)) *((LFLOAT) a(k, j));
  452.           cc(i, j) = (FLOAT) sum;
  453.    }
  454.    *this = cc;
  455.    return *this;
  456. }
  457.  
  458. Matrix &Matrix::operator *= (FLOAT b)
  459. {
  460. #ifndef NO_CHECKING
  461.    Garbage("Matrix *=: *this is garbage");
  462. #endif
  463.  
  464.    for (int i = 1; i <= r; i++)
  465.       for (int j = 1; j <= c; j++) (*this)(i, j) *=  b;
  466.  
  467.    return *this;
  468. }
  469.  
  470. //////////////// elementwise multiplication
  471. Matrix operator % (Matrix &a, Matrix &b)
  472. {
  473. #ifndef NO_CHECKING
  474.    a.Garbage("Matrix %: a is garbage");
  475.    b.Garbage("Matrix %: b is garbage");
  476. #endif
  477.  
  478.    if (a.r != b.r || a.c != b.c)
  479.       a.Nrerror("matrices do not conform in elementary multiplication");
  480.  
  481.    Matrix c(a.r, a.c);
  482.    for (int i = 1; i <= c.r; i++)
  483.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) *b(i, j);
  484.    c.release();
  485.    return c;
  486. }
  487.  
  488. Matrix operator % (Matrix &a, FLOAT b)
  489. {
  490. #ifndef NO_CHECKING
  491.    a.Garbage("Matrix %: a is garbage");
  492. #endif
  493.  
  494.    Matrix c(a.r, a.c);
  495.    for (int i = 1; i <= c.r; i++)
  496.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) *b;
  497.    c.release();
  498.    return c;
  499. }
  500.  
  501. Matrix operator % (FLOAT b, Matrix &a)
  502. {
  503. #ifndef NO_CHECKING
  504.    a.Garbage("Matrix %: a is garbage");
  505. #endif
  506.  
  507.    Matrix c(a.r, a.c);
  508.    for (int i = 1; i <= c.r; i++)
  509.       for (int j = 1; j <= c.c; j++) c(i, j) = a(i, j) *b;
  510.    c.release();
  511.    return c;
  512. }
  513.  
  514. Matrix &Matrix::operator %= (Matrix &a)
  515. {
  516. #ifndef NO_CHECKING
  517.    a.Garbage("Matrix %=: a is garbage");
  518.    Garbage("Matrix %=: *this is garbage");
  519. #endif
  520.    if (a.r != r || a.c != c)
  521.       a.Nrerror("matrices do not conform in elementwise multiplication");
  522.  
  523.    for (int i = 1; i <= r; i++)
  524.       for (int j = 1; j <= c; j++) (*this)(i, j) *=  a(i,j);
  525.  
  526.    return *this;
  527. }
  528.  
  529. Matrix &Matrix::operator %= (FLOAT b)
  530. {
  531. #ifndef NO_CHECKING
  532.    Garbage("Matrix %=: *this is garbage");
  533. #endif
  534.  
  535.    for (int i = 1; i <= r; i++)
  536.       for (int j = 1; j <= c; j++) (*this)(i, j) *=  b;
  537.  
  538.    return *this;
  539. }
  540.  
  541. ///// elementwise division
  542. Matrix operator / (Matrix &a, Matrix &b)
  543. {
  544. #ifndef NO_CHECKING
  545.    a.Garbage("Matrix /: a is garbage");
  546.    b.Garbage("Matrix /: b is garbage");
  547. #endif
  548.  
  549.    if (a.r != b.r || a.c != b.c)
  550.       a.Nrerror("matrices do not conform in division");
  551.  
  552.    Matrix c(a.r, a.c);
  553.    for (int i = 1; i <= c.r; i++)
  554.       for (int j = 1; j <= c.c; j++)
  555.          if (b(i, j)) c(i, j) = a(i, j) / b(i, j);
  556.          else b.Nrerror("zero divisor in elementwise division");
  557.    c.release();
  558.    return c;
  559. }
  560.  
  561. Matrix operator / (Matrix &a, FLOAT b)
  562. {
  563. #ifndef NO_CHECKING
  564.    a.Garbage("Matrix /: a is garbage");
  565. #endif
  566.  
  567.    if (b == 0.0) a.Nrerror("Matrix /: b is zero");
  568.  
  569.    Matrix c(a.r, a.c);
  570.    for (int i = 1; i <= c.r; i++)
  571.       for (int j = 1; j <= c.c; j++)
  572.          c(i, j) = a(i, j) / b;
  573.    c.release();
  574.    return c;
  575. }
  576.  
  577. Matrix operator / (FLOAT b, Matrix &a)
  578. {
  579. #ifndef NO_CHECKING
  580.    a.Garbage("Matrix /: a is garbage");
  581. #endif
  582.  
  583.    Matrix c(a.r, a.c);
  584.    for (int i = 1; i <= c.r; i++)
  585.       for (int j = 1; j <= c.c; j++)
  586.          if (a(i, j)) c(i, j) = b / a(i, j);
  587.          else
  588.             a.Nrerror("zero divisor in elementwise division");
  589.    c.release();
  590.    return c;
  591. }
  592.  
  593. Matrix &Matrix::operator /= (Matrix &a)
  594. {
  595. #ifndef NO_CHECKING
  596.    a.Garbage("Matrix /=: a is garbage");
  597.    Garbage("Matrix /=: *this is garbage");
  598. #endif
  599.    if (a.r != r || a.c != c)
  600.       a.Nrerror("matrices do not conform in elementwise division");
  601.  
  602.    for (int i = 1; i <= r; i++)
  603.       for (int j = 1; j <= c; j++)
  604.          (*this)(i, j) /= ( a(i,j) ) ? a(i,j) :
  605.                     Nrerror("zero divisor in elementwise division");
  606.  
  607.    return *this;
  608. }
  609.  
  610. Matrix &Matrix::operator /= (FLOAT b)
  611. {
  612. #ifndef NO_CHECKING
  613.    Garbage("Matrix /=: *this is garbage");
  614. #endif
  615.    if (!b) Nrerror("zero divisor in elementwise division");
  616.  
  617.    for (int i = 1; i <= r; i++)
  618.       for (int j = 1; j <= c; j++) (*this)(i, j) /=  b;
  619.  
  620.    return *this;
  621. }
  622.  
  623. /////////////////// matrix functions
  624.  
  625. //// sweep out rows k1 through k2,  operates on "a"
  626. Matrix Sweep(int k1, int k2, Matrix& a)
  627. {
  628. #ifndef NO_CHECKING
  629.    a.Garbage("Sweep: a is garbage");
  630. #endif
  631.  
  632.    LFLOAT eps = 1.0e-8, d;
  633.    int i, j, k, n, it;
  634.  
  635.    if (a.c != a.r)
  636.       a.Nrerror("sweep: matrix a not square");
  637.  
  638.    n = a.r;
  639.    if (k2 < k1) { k = k1; k1 = k2; k2 = k; }
  640.    
  641.    for (k = k1; k <= k2; k++) {
  642.       if (fabs(a(k, k)) < eps)
  643.          for (it = 1; it <= n; it++) 
  644.              a(it, k) = a(k, it) = 0;
  645.       else {
  646.          d = 1.0 / a(k, k);
  647.          a(k, k) = d;
  648.          for (i = 1; i <= n; i++) { if (i != k)
  649.             a(i, k) = a(i, k) *(LFLOAT) - d; }
  650.          for (j = 1; j <= n; j++) { if (j != k)
  651.             a(k, j) = a(k, j) *(LFLOAT) d; }
  652.          for (i = 1; i <= n; i++) {
  653.             if (i != k) {
  654.                for (j = 1; j <= n; j++) {
  655.                   if (j != k)
  656.                      a(i, j) = a(i, j) + a(i, k) *a(k, j) / d;
  657.                } // end for j
  658.             } // end for i != k 
  659.          } // end for i
  660.       } // end else
  661.    } // end for k
  662.    return a;
  663. }   /*sweep*/
  664.  
  665. ///// inversion
  666. Matrix Inv(Matrix &a)
  667. {
  668. #ifndef NO_CHECKING
  669.    a.Garbage("Inv: a is garbage");
  670. #endif
  671.    
  672.    if (a.c != a.r)
  673.       a.Nrerror("INVERSE: matrix a not square");
  674.    
  675.    Matrix b(a.r, a.c);
  676.    b = a;
  677.    Sweep(1, b.r, b);
  678.    b.release();
  679.    return b;
  680.  
  681. }   /* inverse */
  682.  
  683.  
  684. ///////////////////////// functions
  685.  
  686. /////// transpose
  687. Matrix Tran(Matrix &a)
  688. {
  689. #ifndef NO_CHECKING
  690.    a.Garbage("Tran: a is garbage");
  691. #endif
  692.    
  693.    Matrix c(a.c, a.r);
  694.    for (int i = 1; i <= c.r; i++)
  695.       for (int j = 1; j <= c.c; j++) c(i, j) = a(j, i);
  696.    c.release();
  697.    return c;
  698. }
  699.  
  700.  
  701. ////// fill with a constant
  702. Matrix Fill(int r, int cc, FLOAT d)
  703. {
  704.    if (r < 1 || cc < 1) {
  705.       Matrix exitnow;
  706.       exitnow.Nrerror("invalid indices in Fill");
  707.    }
  708.    
  709.    Matrix c(r, cc);
  710.    for (int i = 1; i <= r; i++)
  711.       for (int j = 1; j <= cc; j++) c(i, j) = d;
  712.    c.release();
  713.    return c;
  714. }
  715.  
  716. //////// identity
  717. Matrix Ident(int n)
  718. {
  719.    if (n < 1) {
  720.       Matrix exitnow;
  721.       exitnow.Nrerror("n < 1 in identity matrix");
  722.    }
  723.    Matrix c(n, n);
  724.    for (int i = 1; i <= c.r; i++)
  725.       for (int j = 1; j <= c.c; j++) c(i, j) = (i == j) ? 1 : 0;
  726.    c.release();
  727.    return c;
  728. }
  729.  
  730. /////// extract a submatrix
  731. Matrix Submat(Matrix &a, int lr, int lc, int r, int c)
  732. {
  733. #ifndef NO_CHECKING
  734.    a.Garbage("Submat: a is garbage");
  735. #endif
  736.    if ((r > a.r) || (c > a.c)) {
  737.       Matrix exitnow;
  738.       exitnow.Nrerror("SUBMAT: index out of range");
  739.    }
  740.    int r2 = r - lr + 1;
  741.    int c2 = c - lc + 1;
  742.    Matrix b( r2, c2);
  743.    int lrm1 = lr - 1;
  744.    int lcm1 = lc - 1;
  745.    for (int i = 1; i <= r2; i++) {
  746.       for (int j = 1; j <= c2; j++) b(i, j) = a(lrm1 + i, lcm1 + j);
  747.    }
  748.    b.release();
  749.    return b;
  750. }   /* submat */
  751.  
  752. /////  horizontal concatination
  753. Matrix Ch(Matrix &a, Matrix &b)
  754. {
  755. #ifndef NO_CHECKING
  756.    a.Garbage("Ch");
  757.    b.Garbage("Ch");
  758. #endif
  759.    
  760.    if (a.r != b.r)
  761.       a.Nrerror("CH: matrices have different number of rows");
  762.    
  763.    Matrix c(a.r, a.c+b.c);
  764.    
  765.    for (int i = 1; i <= a.r; i++) {
  766.       for (int j = 1; j <= a.c; j++)
  767.          c(i, j) = a(i, j);
  768.    }
  769.    int ii;
  770.    for (i = 1, ii = 1; i <= b.r; i++, ii++) {
  771.       for (int j = 1, jj = 1; j <= b.c; j++, jj++)
  772.          c(ii, jj + a.c) = b(i, j);
  773.    }
  774.    c.release();
  775.    return c;
  776. } /* ch */
  777.  
  778. ///// vertical concatination
  779. Matrix Cv(Matrix &a, Matrix &b)
  780. {
  781. #ifndef NO_CHECKING
  782.    a.Garbage("Cv");
  783.    b.Garbage("Cv");
  784. #endif
  785.    
  786.    if (a.c != b.c)
  787.       a.Nrerror("CV: matrices have different number of columns");
  788.    
  789.    Matrix c( a.r+b.r, a.c);
  790.    
  791.    for (int i = 1; i <= a.r; i++) {
  792.       for (int j = 1; j <= a.c; j++)
  793.          c(i, j) = a(i, j);
  794.    }
  795.    int ii;
  796.    for (i = 1, ii = 1; i <= b.r; i++, ii++) {
  797.       for (int j = 1, jj = 1; j <= b.c; j++, jj++)
  798.          c(ii + a.r, jj) = b(i, j);
  799.    }
  800.    c.release();
  801.    return c;
  802. } /* cv */
  803.  
  804. /////////////// io
  805.  
  806. //////// read a matrix
  807. Matrix Reada(char *filename)
  808. {
  809.    ifstream fp(filename);
  810.    int i, j, r, c;
  811.    char s[241];
  812.    FLOAT f;
  813.  
  814.    if (!fp) {
  815.       cerr << "Cannot open " << filename << ".\n";
  816.       Matrix exitnow;
  817.       exitnow.Nrerror("ReadA: file open failure");
  818.    }
  819.    fp.getline(s, 240, '\n');
  820.    fp >> r >> c;
  821.    Matrix a(r, c);
  822.    for (i = 1; i <= r; i++) {
  823.       for (j = 1; j <= c; j++)
  824.          if (fp >> f) a(i, j) = f;
  825.          else a.Nrerror("Reada: error reading input");
  826.    }
  827.    fp.close();
  828.    a.release();
  829.    return a;
  830. }
  831.  
  832. /////////// write a matrix
  833. void Writea(char *filename, Matrix &a, const char *filetitle)
  834. {
  835. #ifndef NO_CHECKING
  836.    a.Garbage("Writea: a is garbage");
  837. #endif
  838.    ofstream fp(filename);
  839.    int i, j;
  840.  
  841.    if (!fp) {
  842.       cerr << "Cannot open " << filename << ".\n";
  843.       a.Nrerror("Writea: file read failure");
  844.    }
  845.    fp << filetitle << "\n";
  846.    fp << a.r << " " << a.c << "\n";
  847.    for (i = 1; i <= a.r; i++) {
  848.       for (j = 1; j <= a.c; j++)
  849.          if (!(fp << a(i, j) << " "))
  850.             a.Nrerror("error writing output");
  851.          fp << "\n";
  852.    }
  853.    fp.close();
  854. }
  855.  
  856.